Dense granular flow around a penetrating object: 
Experiments and hydrodynamic model 
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We present in this Letter experimental results on the bidimensional flow field around a cylinder 
penetrating into dense granular matter together with drag force measurements. A hydrodynamic 
model based on extended kinetic theory for dense granular flow reproduces well the flow localization 
close to the cylinder and the corresponding scalings of the drag force, which is found to not depend 
on velocity, but linearly on the pressure and on the cylinder diameter and weakly on the grain size. 
Such a regime is found to be valid at a low enough "granular" Reynolds number. 

PACS numbers: 45.70.-n, 45.50.-j 



Introduction - Describing the motion of an obstacle 
through granular material is the subject of recent and 
intensive research with applications to industrial config- 
urations but also to biological and earth science with, ex- 
ample include animal locomotion in sand [l| and impact 
cratering Q . If the motion of an object in a simple fluid 
is known for a long time, especially in the viscous regime 
at low Reynolds number where the fluid flow and the drag 
force are analytically known since Stokes' calculation, the 
motion of an object in granular matter is still an open 
question. Such a problem if of fundamental interest, with 
numerous open questions of statistical physics concerning 
(for instance) the solid-liquid or jamming transition Q. 
Numerous studies have been done concerning the drag 
force on an object in vertical or horizontal motion in gran- 
ular matter [441 lj. All these studies find that the drag 




force does not depend on the velocity at low velocities, is 
proportional to the size of the object, and to its depth. 
As in hydrodynamics, the drag force has been shown to 
depend on the exact shape of the object ||, and also 
vertical lift forces can develop during horizontal motion 
[l2T ] . Flow observations of grains have been also reported 
in chute flow around a fixed cylinder 13|, [L4| and in the 
two-dimensional situation of a disk pulled at a constant 
force in a horizontal assembly of disks on a vibrated plate 
Fluctuations have been observed in the force or in 
the velocity with some "stick-slip" behavior in some cases 
d, 0] , and the force may depend crucially on the pack- 
ing volume fraction [3|— 15[| - In this Letter, we investigate 
by Particle Image Velocimetry (PIV) measurements the 
flow around a cylinder penetrating into a dense granular 
packing together with force measurements. By a contin- 
uum hydrodynamic model based on the kinetic theory 
extended to dense granular systems, we recover the ex- 
perimental results of the shear localization close to the 
cylinder with the view of a "hot" cylinder in motion in a 
viscosity-dependent temperature fluid, together with the 
good scaling for the drag force. 

Experiments. -The experiments consist in a horizontal 



FIG. 1. Typical instantaneous velocity field obtained by PIV 
measurements for a cylinder of diameter d — 20 mm pene- 
trating in glass beads of diameter d g — 1 mm, at the velocity 
Vo = 5 mm s _1 at the depth Zb = 65 mm. 



steel cylinder of diameter 10mm< d <40mm and length 
b = 40 mm penetrating in a rectangular box (0.1m in 
height, 0.2 m in width and 40 mm in thickness) filled 
with monodisperse millimetric glass beads of diameter 
0.5 mm< d g <4mm and density p g — 2.5 10 3 kg m~ 3 . 
The granular medium is prepared by gently stirring the 
grains with a thin rod and the surface is then flatten us- 
ing a straightedge. We have checked that this preparation 
leads to reproducible results with only small variations. 
The solid volume fraction is <£> ~ 0.62 characteristic of a 
dense granular packing and the density of the granular 
medium is thus p — p s <f> ~ 1.5 10 3 kg m~ 3 . The cylin- 
der which is fixed and related to a force transducer by 
a vertical thin rod, is first above the grain surface and 
penetrates gradually into the granular packing as the box 
is raised up by a stepper motor along a vertical transla- 
tion guide at a constant velocity Vo ranging from 0.1 to 
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FIG. 2. (a) Drag force F on the cylinder as a function of 
the penetration depth zt for different cylinder diameters from 
d — 10mm to d — 40mm (d g = 1mm, Vo = 5mm s" 1 ). 
(b) F/zb for the same data as (a), (c) < F/zt > as a function 
of Vo (d = 20 mm, d g = 1 mm) and (d) as a function of d g 
(d — 40mm, Vo = 5mm s _1 ). < F/zt > is the average of 
F/zb over the range d/2 < Zb < 70 mm. 



100 mm s _1 . Very careful alignment is taken to prevent 
any blockage of the cylinder during the motion, and the 
force at the glass walls without grains is totally negligible. 
The front and rear wall of the box are in glass allowing vi- 
sualization of the granular flow around the cylinder. The 
images taken from a fast video camera (up to 1000 im- 
ages per second in the full resolution 1024 x 1024 pixels) 
are analyzed by a PIV software to get the velocity field 
of the grains. As the camera is fixed in the laboratory 
frame together with the cylinder, the obtained velocity 
field shown in Fig. [T] is the velocity of the grains in the 
frame of reference of the cylinder. 

The measured drag force on the cylinder is observed to 
increase with the depth Zh during its penetration (Fig. 2a) 
with a ratio F/zb constant to ±10% over the range 
d/2 < zi, < 70mm (Fig. 2b). The force is found pro- 
portional to the cylinder diameter (Fig. 2b) and roughly 
independent of the velocity (Fig. 2c). We find also a non- 
linear dependence of the force with the grain diameter: 
The force is about constant at large enough grain size 
iAg ^ 1 mm) but increases with decreasing grain size 
(d g < lmm)(Fig.2d). 

Concerning now the grain flow field around the totally 
immersed cylinder, the first important and key obser- 
vation is that the velocity field is stationary during the 
penetration process, meaning that it depends neither on 
the depth nor on the increasing granular pressure. The 
velocity field v(x, z) can thus be averaged during each 
penetration run to extract at each point the mean local 
velocity v(x, z) whose time fluctuations can be related to 
the local granular temperature T, by T =< (v — v) 2 >. 
Due to the geometrical configuration, we have used the 
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FIG. 3. Radial velocity functions A r and Ag as a function of 
the distance r from the cylinder. Experimental data (o) A r 
and (A) Ag from PIV measurements; ( — ) Profiles obtained 

by using granular kinetic theory; ( ) Classical Newtonian 

profiles. 

most appropriate cylindrical coordinates (r, 9) where r is 
the radial distance from the cylinder center and 9 is the 
angle relative to the downward z-axis of motion (with 
thus 9 = down): v(r, 6) = v r (r,9)e r + vg(r,9)eg, with 
the radial and azimuthal components of the velocity v r 
and vg. As in the classical hydrodynamics situation of 
a Newtonian fluid, we have checked that v r and vg can 
be decomposed into cosine and sine functions of 9 and 
radial functions A r (r) and Ag(r) : v r = — VoA r (r) cos 9 
and vg = VoAg(r) sin 9. The radial functions A r (r) and 
Ag(r) extracted from measurements in the azimuthal 
range — tt/2 < 9 < ir/2 are displayed in Fig. [3] and show 
a strong shear localization when compared to the classi- 
cal viscous Newtonian case, with exponential variations 
scaled by the cylinder diameter. At a distance from the 
cylinder surface larger than about one cylinder diameter 
(r > 40 mm in Fig. [3]), the grain velocity vanishes. The 
overshoot of Ag(r) expected from mass conservation in 
the bidimensional configuration is localized close to the 
cylinder and relaxes to the asymptotic value 1 (no grain 
flow) with an inflexion point in the present case in con- 
trast with the long-range decay of the Newtonian case at 
Re = (Couette-Poiseuille form). 

A typical radial profile of the granular temperature 
T(r) along the vertical downward line (6 = 0) is displayed 
in Fig. [4j A domain of roughly constant temperature can 
be seen around the cylinder with the plateau value To 
followed by an exponential decrease. The temperature 
profile is roughly the same for different values of 9 and 
the plateau value To is found to vary with the penetration 
velocity as T ~ V 2 ( see mset °f Fig. EH- 

Hydrodynamic model - The granular mean velocity and 
its fluctuation show important spatial variations. The 
strong localization of the granular flow that we observe 
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FIG. 4. (o) Experimental and ( — ) theoretical profiles of the 
granular temperature T/Vq as a function of the distance r 
from the cylinder, for 8 — 0. Inset: (o) experimental and 
(■) theoretical plateau temperature To as a function of Vb; 
( ) power fit of the data of the form To oc Vq. 



is a rather usual feature of disordered matter where shear 
bands are commonly observed (see. (TH[| for a recent re- 
view). We show here that the observed shear bands may 
be understood using an hydrodynamic description. The 
starting points are that local momentum and energy bal- 
ance equations, written as 161 ]: 



dT 

p-rr = c : k — V • q — eT, 



(1) 
(2) 



where p is the effective fluid density, a the stress ten- 
sor, q the heat flux, k the velocity-gradient tensor, and 
e the temperature-loss coefficient. The stress and the 
heat flux are related to the velocity, temperature and 
pressure by phenomenological equations. The choice of 
these phenomenological equations to describe granular 
matter is still matter of debate. If momentum transfer 
and dissipation occur during binary collisions between 
grains, granular material may be treated using an inelas- 
tic gas theory 17[, leading to a Newtonian fluid with 
(7 = —PI + 2rjK, and Fourier's law q = — AVT for the 
heat flux. Here P is the pressure, r\ the viscosity, A the 
thermal conductivity and I the unit tensor. For simplic- 
ity, we neglect dissipation and heat transport associated 
with compressibility and we treat granular matter as an 
incompressible fluid as experimentally observed. A 'pri- 
ori, in a dense granular flow, non-binary collisions cannot 
be neglected and transport coefficients from inelastic gas 
theory Tjj are no longer valid. By analogy with glassy 
materials, modification of the viscosity divergence near 
jamming has been suggested (H)]. Numerical simulations 
of 2D granular materials seem then to show viscosity di- 
vergence at packing fractions lower than random close 



packing [l6[ [ijj]. Since in our experiment, most of the 
shear is located in a region of high velocity fluctuations, 
we are not very close to random close packing, and there- 
fore we use the Enskog expressions of the phenomenolog- 
ical coefficients. We obtain tl7l : 



P = Pg Tf P {$) 
\ = p g d g VTf x ($) 



e = (i-0 



(3) 
(4) 
(5) 

(6) 



where fp, f v , f\ and f s are non-dimensional functions 
of the solid volume fraction $, and e is the velocity 
restitution coefficient. For <f> > 0.5, those functions 
vary with the same dependence on $ 17J. So we have 
r,(P,T) ~ rj x {Pd g /VT), X(P,T) ~ A x (Pd g /VT) 
and e(P, T) ~ e x (P/d g y/T), with i l0 ~ 0.28, A ~ 1.06 
and £q ~ 0.34 [13] with a standard value e ~ 0.9 for 
glass beads. Those expressions of transport and dissipa- 
tion coefficients emphasize the fact that 4> is not fixed 
in our experiment but may vary slightly from point to 
point in response to pressure and granular temperature 
variations. 

The momentum and heat equations (11I2[) with T — P 
dependent transport coefficients are solved numerically 
for the stationary flow around a cylinder located at the 
center of a L x L square box. The momentum equation 
is solved using a Lattice-Bolzmann solver (BGK based 
D2Q9 model) with non-slip velocity conditions on the 
cylinder, constant pressure Pq at the top of the square 
box, and constant upward velocity Vo on the other 
sides. The heat equation is solved using a finite difference 
scheme with the condition e r ■ VT = on the cylinder 
and T = on the sides of the square box in agreement 
with experimental findings. The transport coefficients 
are taken initially homogeneous in the box and the veloc- 
ity field is first computed by solving momentum equation 
with these initial values. From the obtained pressure and 
velocity fields, the source of heat a : n is calculated and 
the heat equation is then computed leading to a new tem- 
perature field T(x, z). With the new corresponding fields 
of transport coefficients as inputs, the momentum equa- 
tion and then heat equation are solved again and so on. 
After a few iterations, stationary temperature, pressure 
and velocity fields which verify (jl!2p are obtained. In 
order to prevent numerical instabilities, viscosity is kept 
in a finite range. We have checked that the stationary 
solution is not sensitive to the initial guess of temper- 
ature, neither to the cut-off of viscosity. We have also 
checked that other boundary conditions at the cylinder 
(partial velocity slip and Robin condition dT/dr oc T for 
the temperature) do not change significantly the velocity 
and temperature fields. 

Results— Our hydrodynamic model reproduces quite 
well the experimental features observed in the experi- 
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merits. The radial temperature profile shows the same 
shape as in the experiments, with a plateau of high tem- 
perature T close to the cylinder followed by an expo- 
nential decrease (Fig. 0J. The temperature plateau is 
also found to be proportional to Vq as in the experi- 
ments (see inset of Fig. 4) . The computed velocity field 
is found close to the experimental one as shown in Fig. [3] 
with a shear localization near the cylinder, and is far 
from the classical Newtonian case. This shear localiza- 
tion may then be interpreted as a consequence of the 
strong coupling between viscosity and temperature. At 
a given pressure the viscosity varies as rj ~ 1/VT. So, at 
a given viscous shear stress a v , the production of heat is 
proportional to ~ y/T '. This mechanism creates a 

self-lubricating layer of low viscosity near the cylinder. 

The drag force on the cylinder is calculated by inte- 
grated the stresses on the disk, taking into account both 
the "pressure" term (from normal stresses) and "viscos- 
ity" term (from shear stresses). If these two terms are 
equal in the Newtonian case, the pressure term is here 
about twice the viscosity term. The total calculated drag 
force is found to be independent of the velocity, propor- 
tional to the cylinder diameter and to the pressure, in 
agreement with the experimental observations by con- 
sidering that pressure is proportional to depth. That 
dependance may be understood in the hydrodynamical 
model by considering a "granular" Reynolds number Re 
= pVod/rj w based on the viscosity r] w near the cylinder. 
We have checked that all these previous findings corre- 
spond to a low Reynolds number regime (Re < 1). In low 
Re hydrodynamics, one expects that the force scales as 
i]w x Vq, and then varies here linearly with the pressure, 
and independently of the velocity as r\ w is proportional to 

1/2 1 

the pressure and to T Q ' thus to V . We also found 
numerically the same non linear variation of the drag 
force with the grain size as in the experiments. Note 
that this variation is hard to infer simply from the set 
of equations. When Re > 1, the velocity field no longer 
exhibits up/downstream symmetry, and the pressure and 
temperature profiles are also different from the Re < 1 
case. 

Concluding remarks.— We have investigated experimen- 
tally the penetration of a cylinder at a constant veloc- 
ity inside a dense granular packing with both force and 
velocity field measurements, and we have modeled this 
problem by a continuum hydrodynamic approach. The 
finding of a strong shear localization close to the cylinder 
can be viewed by the coupling between viscosity and tem- 
perature in the problem of a self-heated cylinder. The 
localization of the flow near a sedimenting hot sphere 
has indeed already been reported in classical fluids with 
temperature dependant viscosity [2(j. Such a shear lo- 
calization has been seen also for a sphere sedimenting in 
a non Newtonian fluid with shear thinning behavior [21| . 
The experimental findings of a force regime independent 
of velocity, proportional to the depth and to the cylinder 



diameter are recovered by our model based on kinetic the- 
ory adapted for dense granular systems, and have been 
shown to correspond to an hydrodynamic regime of low 
"granular" Reynolds number. Other models exist for 
dense granular flows such as the one based on a local 
rheology with a friction coefficient depending on a 
non-dimensional shear rate / [22} . Such models may also 
lead to the observed flow localisation since it corresponds 
to visco-plastic/shear thinning behavior. Another issue 
to explore would be the different force value measured 
in the plunging and the withdrawal situations 0, lioj . 
with the role played by gravity and the boundary con- 
ditions (bottom wall vs free surface). In addition, the 
force scaling is expected to change in a higher Reynolds 
regime from a "viscous" scaling to an "inertial" scaling 
which may explain the complex force terms measured 
in impact situations @, [H , 24 1 where the "granular hy- 
drodynamic" regime change certainly from high to low 
Re during the penetration process. Applications to non- 
stationary granular flows may also been considered. 
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